Infiltration through porous media 
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We study the kinetics of infiltration in which contaminant particles, which are suspended in a 
flowing carrier fluid, penetrate a porous medium. The progress of the "invader" particles is impeded 
by their trapping on active "defender" sites which are on the surfaces of the medium. As the 
defenders are used up, the invader penetrates further and ultimately breaks through. We study this 
process in the regime where the particles are much smaller than the pores so that the permeability 
change due to trapping is negligible. We develop a family of microscopic models of increasing 
realism to determine the propagation velocity of the invasion front, as well as the shapes of the 
invader and defender profiles. The predictions of our model agree qualitatively with experimental 
results on breakthrough times and the time dependence of the invader concentration at the output. 
Our results also provide practical guidelines for improving the design of deep bed filters in which 
infiltration is the primary separation mechanism. 



I. INTRODUCTION 

In depth filtration, suspended particles in a fluid are 
removed during their passage through a porous medium 
|H,p|. The basic dynamics of depth filtration is deter- 
mined primarily by the pore structure of the filter, the 
particle size distribution, and by various physicochemical 
and hydrodynamic details. If the particle size is larger 
than the typical pore size, particles get stuck relatively 
quickly. The permeability of the filter decreases steadily 
during this process and drops to zero when clogging is 
reached. This process is often referred to as sieving, or 
straining |3|. Conversely, if particles are much smaller 
than the pore size and if particles are trapped only at 
the interfaces of the porous medium, the flow field is only 
slightly affected by the trapping. The goal of this paper is 
to provide a general understanding of this latter process 
of infiltration by microscopic network modeling. 
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FIG. 1. Idealized picture of infiltration. Suspended parti- 
cles are trapped at "defender" sites on pore surfaces. Once 
defenders are occupied, subsequent particles pass by freely. 

Infiltration underlies many practical situations, such 
as underground waste disposal B, gas mask design, or 
drinking water filters [g . Typically, sub-micron size con- 
taminant particles are suspended in a carrier fluid and 
flow through a porous material, such as a sand filter 
whose typical grain size is much larger than the con- 
taminant particles, or an ion exchange filter ||5|] where 
the contaminant size is molecular in scale. In such cases, 
one can neglect the change of the flow field due to par- 
ticle trapping B, an approximation which considerably 
simplifies theoretical analysis. 



The kinetics of infiltration is controlled by the micro- 
scopic mechanisms for the trapping of the invader par- 
ticles. Typically each pore can hold a limited number 
of particles due to a finite surface area or a finite range 
of the surface potential. When all the available surface 
area is covered by particles, subsequent invaders flow pas- 
sively through the filter without being trapped. Our ba- 
sic goal is to understand the kinetics of this infiltration 
and the ultimate breakthrough of the invader, as well as 
the evolution of the invader and defender density profiles 
as functions of downstream position and time. 

Previous work on infiltration in porous media has often 
been based on a macroscopic convection-diffusion equa- 
tion description, with reaction terms introduced to ac- 
count for particle trapping |H,pln|. Another approach 
has been to use a single absorbing sphere to calculate 
the collection efficiency at the initial stage of filtration 
[^. While numerical simulations of these models have 
some predictive power, it is hard to develop a connection 
between this macroscopic approach and basic features of 
the microscopic process, such as the concentration pro- 
files of the trapped and fiowing particles. 

For filtration by straining, models based on a discrete 
network description of the filter medium are relatively 
well developed ||9|-p^. To our knowledge, however, there 
has been no microscopic network modeling work on infil- 
tration. As in the case of straining, a spatial density gra- 
dient naturally arises in infiltration, since particles begin 
to deposit at the upstream end of the filter and advance 
downstream as the filter gets used up. The density gra- 
dient is experimentally observed as the time-dependent 
output concentration g]. In this paper, we will account 
for this basic experimental observation by using a dis- 
crete network approach. 

Practical questions raised by infiltration are the break- 
through time, which is defined as the time for the output 
concentration to reach a specified threshold level, and the 
filter efficiency, which is related to the fraction of the fil- 
ter material actually used before breakthrough. Clearly, 
it is desirable to use as much of the filter material as 
possible before breakthrough occurs. 



This paper is organized as follows. In Sec. ||, we intro- 
duce the basic parameters that govern particle trapping 
and provide a qualitative picture of infiltration. In the 
following sections, we construct a sequence of discrete 
models with increasing complexity and realism to ulti- 



mately provide a lattice network description. In Sec. |I|, 
we discuss the case of a one-dimensional (ID) chain of 



trapping sites and in Sec. IV, we analyze infiltration in 
the bubble model to provide a mean-field-like descrip- 
tion. Building on these results, we then turn to simu- 
lations of infiltration on tube lattice networks in Sec. M. 
We summarize and compare our results with experiments 
in Sec. Ivl. 



II. BASIC PICTURE 

The two basic characteristics of particle trapping are 
the efficiency of an unoccupied trapping site and the 
number of trapping sites in a pore. We introduce the 
trapping probability 7 as the probability that a particle 
is trapped upon encountering an open collector site. The 
parameter 7 thus represents the strength of the particle- 
collector interaction and accounts for the possibility that 
contact between particles and the filter grains may not 
necessarily lead to deposition ||lj]. While this simplifies 
the complicated adsorption mechanism, later we show 
that the basic feature such as the invasion front propa- 
gation velocity (Fig.0) is independent of the interaction 
details. 

Next we introduce the capacity Cx{t) as the number 
of particles a pore at position x can hold at time t. In 
the case of non-coagulating particles which cannot get 
trapped on top of an already adsorbed particle, the ini- 
tial capacity is proportional to the inner surface area of 
the pore and then decreases as the pore surface is covered 
by particles. For simplicity, we ignore multiple trapping 
on an already occupied collector site as well as particle 
re-launching. The key factors which determine the dy- 
namic behavior of the system are geometric, such as the 
capacity of a clean filter and the pore size distribution, 
and kinematic, such as the particle concentration and the 
fiow rate. More refined models for particle trapping can 
be incorporated within our basic modeling. 

Consider a generic infiltration process based on the 
above concepts. Initially a layer in the clean filter has 
a total capacity Ctot independent of the downstream po- 
sition. At t = 0, a fluid which contains a mixture of 
invader and non-reacting tracer particles enters the filter 
whose flow rate is determined by the steady-state solu- 
tion of d'Arcy's law. Tracer particles passively follow the 
fluid motion and advance with the average flow velocity 
vq. The width of the tracer density proflle spreads as t^/^ 
due to hydrodynamic dispersion (Fig. 0). Invader parti- 
cles flrst encounter clean collector sites. Because each 
such encounter leads to deposition with probability 7, 
the survival probability of the particles in this leading 



invaded region decreases exponentially with downstream 
position as illustrated in Fig. || As particles advance 
and get trapped, the pore capacity decreases and subse- 
quent particles are more likely to survive, giving rise to 
an advancing invasion front with a velocity v < vo- In 
principle, the propagation velocity and shape of the front 
are functions of time. However, at long times these fea- 
tures approach steady-state values. Another important 
feature is that the trailing edge of the capacity proflle 
decays as a power of the distance |^| (^ < 0) from the 
invasion front whose location is deflned, e.g., as the po- 
sition where c = Ctot/2. For large |^|, any reasonable 
deflnition for the front location can be used. 

The existence of different propagation velocities, vq for 
the pure fluid and v < vq for the contaminant, leads to 
purification of the liquid. The filter can be used until 
the invaded region reaches the outlet end. For a filter of 
length L, the breakthrough time will be of the order of 
L/v, so the amount of throughput will be approximately 
proportional to Lvq/v. 
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FIG. 2. Infiltration profiles. Horizontal direction is down- 
stream. The invader density profile (shaded) is exponential 
in the invaded region, while the capacity profile has a power 
law tail, c~ l^l""^. 



III. ONE-DIMENSIONAL MODEL 

As a preliminary, we study infiltration in a one- 
dimensional chain of identical pores at k = 0,1,2,.... 
First we consider the case where each pore can accommo- 
date only one particle and then we generalize to multiple 
capacity pores. 



A. Single capacity pores 

We choose a time unit such that one particle is in- 
jected at each discrete time step. Multiple particle injec- 
tion leads to a different particle density and will affect 
only the overall scale factor and not change qualitative 
features of the system. The carrier fluid advances by 
one pore distance at each time step; that is, its veloc- 
ity is unity. When particle trapping occurs in a pore 
at fc, the capacity Ck{t) changes permanently from 1 to 
0. At time i, a particle at pore k gets trapped in that 
pore with probability 7 if Ck(t) = 1. A particle advances 
to the next pore in one time step with probability 1 if 



Ck{t) = 0. Based on these elemental steps, we introduce 
the following two probability densities: 

• pk (t) : The probability that a freely- moving particle 
is in pore k at time t. 

• Qkit)'- The probability that the pore at site k is 
unoccupied; that is Cfc(i) = 1. 

The corresponding master equations for Pk{t) and qk{t) 



Pk{t + 1) =pk-i{l -79fc_i), 
qk{t + 1) = gfe(l-7Pfc), 



(1) 
(2) 



where we drop the argument t on the right hand side for 
simplicity. Unless there is a possibility for confusion, we 
will not write the argument t in related formulae. Since 
a particle advances to the next pore in one time step, 
Pk{t + 1) depends on pk-i{t). The term (1 — jqk-i) in 
Eq. (Q) is the probability that the particle at fc — 1 does 
not get trapped by an unoccupied pore also at k — 1. 
Similarly, the term (1 — 7Pfc) in Eq. (0) is the probability 
that pore k does not trap a free particle at time t. The 
initial and boundary conditions for these equations are: 

Pfc>i(0) = ft(0) = l; Po{t) = l <Zoo(i)-l (3) 

If the trapping probability 7 is small, a particle can ad- 
vance many pores without being trapped, so that Pk{t) 
and Qk (t) vary slowly in space and time and a continuum 
approximation can be applied. Letting k + 1 —> x -\- Sx 
and t + 1 —^ t + 6t, Eqs. (|l|) and (||) become, to lowest 
order. 



dtp + vodxP 
dtq 



'ipq, 
'jpq, 



(4) 



where vq = Sx/St — 1, and 7 — > j/St is a redefinition of 
the trapping probability in units of the infinitesimal time 
increment. 

In a co-moving reference frame, S^ = x — vt, with v 
the invasion front propagation velocity shown in Fig. S 
(which is yet to be determined), Eqs. (|j) become 



{vo - w)9^p + dtp = -jpq 
- vd^q + dtq = -ipq. 



(5) 
(6) 



Let us first examine the steady state solution of 
Eqs. (g) and (|^). Setting the time derivatives to zero, 
subtracting Eq. (|6|) from Eq. (|5|), and integrating with 
respect to f gives 



{vq - v)p{^) + vq{i) = const. 



(7) 



The integration constant can be determined by applying 
Eqs. (ph in the co- moving frame. As ^ ^ —00, p ^ 1, 
q — > 0, and as ^ — > 00, p — > 0, g — > 1. These immediately 
give V = uo/2 = 0.5. Note that v is determined entirely 
from the boundary condition (and wq) in the co- moving 



frame, and not from the interaction strength 7. This 
feature continues to hold for all the models in this paper. 
Using Eq. (0) with v — wo/2, Eqs. (||) and (||) can be 
solved to give 
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m = 



1 



.-i/io ' 



(8) 



(Fig. H). Thus ^0 = Vo/2^ is the characteristic width w 
of the profile. Notice also that the profiles oip and q are 
symmetric about their intersection. We verified both the 
dependence of the width on 7 and the profile shape pre- 
dicted by Eq. (||) by numerical integration of the master 
equations (M and (0). 
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FIG. 3. Steady state profile in the single capacity model. 

One subtle point is the rate of approach to the steady 
state. First, we find that the asymptotic propagation ve- 
locity V = iiq/S is reached before the asymptotic profile 
is established. This arises because v is determined by 
the boundary conditions, and not by interaction details. 
We also verified this feature numerically. Adopted this 
asymptotic velocity, Eqs. (m and (ra) are then symmetric 
in p and q. In fact, the system is identical to two-species 
annihilation, A + B ^ 0, where each species is ballisti- 
cally injected from opposite sides with velocity +wo/2 for 
the ^s and —vq/2 for the Bs. 

Thus most of the the time variation in Fig. ^ occurs 
in the reactive region of width ly, where w = at < = 0. 
Integrating Eq. (Eh from —^1 <C to ^2 ^ gives. 



Wo 



i2 



pd^^-'f 



«2 



pqd£,. 



«i 



(9) 



The integral on the left hand side is the area under the 
curve p{^) between —^1 and ^2, whose time dependence 
mainly comes from the change in w. On the right hand 
side, the integral is approximately proportional to w, 
since pq is significantly different from zero only in the 
reactive region. We can then rewrite Eq. (0) as 



I'a 



— — + cidtw ~ — 7C2W, 



(10) 



where Ci and C2 are constants. Integrating Eq. (001) and 
applying the condition w{oo) ~ ^ gives an exponential 

decay to the steady state w{t) ^ ^{^ — e~'^^^*). 

It is worth emphasizing that the symmetry between 
the invader and defender is generally responsible for the 



relation v = uo/2. At the inlet, invaders are injected with 
velocity wq, and the invasion front advances with velocity 
V, with one invader particle annihilating with one de- 
fender site. In the reference frame moving with velocity 
vq, the situation is reversed. The invaders are at rest and 
defenders are injected with velocity vq from the opposite 
direction. Therefore the invasion front advances with ve- 
locity Vq — V. Since these two reference frames describe 
the system in the same way, the front velocities should 
be the same; that is, v = vq — v, ot v = vq/2. 



B. Multiple capacity pores 

Now we consider the case where each pore can trap M 
particles, that is, the initial pore capacity is Cfe(O) = M. 
We again follow the previous rules of injecting a single 
particle and advancing a particle by one pore (wq = 1) 
at each time step. Multiple particle injection or differ- 
ent injection intervals again simply changes the overall 
concentration and time scale. 

In a multiple capacity pore, the probability of encoun- 
tering an open trap in a pore needs be considered, in 
addition to the trapping probability 7 upon encounter 
with an open trap. Generally, the encounter probability 
decreases as more particles get trapped, since the inner 
pore surface area available for trapping shrinks. When 
fluid mixing within a pore is weak, a particle can en- 
counter only one trap, either open or occupied. Then 
the encounter probability is approximately proportional 
to the fraction of the open surface area. On the other 
hand, if the mixing is perfect, a particle encounters all 
the available traps in a pore before exiting. 

For practically relevant situations, pores are suffi- 
ciently short so that a particle in a pore follows stream- 
lines without transverse diffusive mixing Pq]. In what 
follows, we consider this limit of weak mixing. For a pore 
with n out of M traps available, the encounter probabil- 
ity is n/M ^ and the overall trapping probability of this 
pore is r„ = ^n/M . In writing this expression, we ignore 
the possibility that a particle far from the pore wall does 
not encounter any traps. In Sec. ^, we argue that this 
volumetric effect does not change the basic behavior of 
infiltration. 

To describe the evolution of the system, we use the 
same single-particle probability density Pkit) as in the 
single-capacity pore system, but modify the probability 
density for the capacity as follows: 

• Qkit): The probability that a pore at position k 
contains n open traps. This is the same as the 
probability that Ck{t) — n, ioi < n < M . 

Following similar reasoning as that applied to deduce 
Eqs. (|l|) and (g), the master equations for pk{t) and 9j!(t) 
are 



M 



Pkit + 1) = Pfc-1 [gLi + Y. gri(l - r„)] 



)i=i 



(11) 



q''kit+l) = qn^-PkTn)+Pkq'^+'T, 



n+l- 



(12) 



In Eq. (jT^), q^_i{t) accounts for the case that the pore 
k — 1 has zero capacity. Other terms in Eq. ( Jll] ) cor- 
respond to cases when the capacity is different from 0, 
with (1 — T„) the survival probability for each case. In 
Eq. (|l^ ) q^{l—pkTn) is the probability that the pore with 
capacity n does not trap a particle, and Pkq^ ^n+i i^ 
the probability that the capacity decreases from n -I- 1 to 
n by a particle trapping event. Hence the last term is 
absent when n = M . 

We simplify Eqs. (^I]) and (|l^) by introducing the av- 
erage capacity of a pore at position k, 



M 



Qk{t)^Y.''ql{t). 



(13) 



n=l 



This gives the average number of sites still available for 
trapping in the pore. Now by multiplying Eq. ( |l2[ ) by n, 

summing from 1 to Af , and using X]n=i 9fe (*) = l~Q'fc(i)i 
we obtain 



Pk{t + I) = Pk-i{l 



M 



^-1) 



h{t + i)^Qk{i-j-^Pk) 



(14) 
(15) 



These are identical in form to Eqs. (Il|) and (0), so the 
same steady state analysis applies. We transform to a 
co-moving frame and take the continuum approximation 
to reduce the rate equations to Eqs. (||) and (g) with 
7 -^ 7/M and q ^ Q- The boundary conditions are also 
the same as in the case of single-capacity pores, except 
Q{^ ^ 00) = M. Combining these results give 



vo 



1 + M 



w(oo) 



VqM 



7(1 + M) ■ 



(16) 



Notice that for M = 1, Eq. ( [Tq ) reduces to the single ca- 
pacity case, while for M — > 00, w — > 0. This means that 
there is no steady state for the case of infinite capacity 
pores. 

We can generalize the symmetry argument given in the 
single capacity case to find the propagation velocity in 
Eq. (|l6|). At the input, the flux of invaders moving with 
the carrier fluid is equal to 1 x wq. Similarly, in the refer- 
ence frame moving with velocity vq , the flux of defenders 
is M X Vq, while the invaders are at rest. Because one 
invader annihilates with one defender, the two particles 
are kinetically indistinguishable. Therefore, if a particle 
flux of 1 X Wo results in a front moving with velocity v, 
the front velocity produced by a flux of M x vq should 
be Mv, which, in turn, equals vq — v in the moving ref- 
erence frame. By this equivalence, Eq. (^^ immediately 
follows. 

In the limit of perfect mixing, a particle encounters all 
traps in the pore. The overall trapping probability with 



n open traps is then T„ = 1 — (1 — 7)". In the hmit of 
smah 7, T„ ~ 771, thus the analysis is exactly the same as 
in the poor mixing case except without the factor \/M in 
T„ . The propagation velocity of the front is the same as 
in Eq. (nq) , since this velocity is independent of trapping 
mechanism, while the width varies as w ^ wo/7(l + M). 
Notice that w is a decreasing function of M . This arises 
because a particle must survive all the traps in a pore 
before advancing to the next pore. Finally, for a mix- 
ing mechanism which is intermediate between the two 
limits of perfect and poor mixing, the propagation ve- 
locity will be vo/(l + -W), while the width of the front 
will lie between the limiting values of ^0/7(1 + M) and 
voM/7(l + M). 



IV. THE BUBBLE MODEL 



particle chooses a tube in the next downstream bubble 
according to flow induced probability <i>(r). 



$(r) = 



(18) 



/dr' /(r') r''' ' 

in which the probability of choosing an outgoing tube of 
radius r is proportional to the flow rate into the tube, r^ 

Since tubes of different flow velocities give the domi- 
nant mechanism for dispersion, the radial dependence of 
the local flow velocity in a tube (Taylor dispersion) is ig- 
nored. Thus we assume that a particle moves with the 
average flow velocity v{r) along the tube. We now inves- 
tigate the hydrodynamic dispersion of passive Brownian 
particles which are carried by the background fluid in the 
bubble model, in the absence of any trapping. This will 
provide the concepts and tools necessary to understand 
infiltration in the bubble model. 



We now study the bubble model as a logical next step 
towards understanding infiltration in porous media. The 
bubble model was introduced to account for the breaking 
of fibers |lq| , extremal voltages in resistor networks [O , 
and later to filtration kinetics [Q. The bubble model 
consists of L "bubbles" in series, each of which is a paral- 
lel bundle of w tubes, with each tube representing a pore 
(Fig. ^ . A bubble can be viewed as a single layer of par- 
allel bonds in a lattice with all the ends "shorted" . This 
model has multiple paths, as in real porous media, and 
is sufficiently simple to be amenable to analytic study. A 
useful feature of the bubble model is that for straining 
dominated filtration, this model predicts similar behavior 
to that of lattice networks llltl . 





FIG. 4. The Bubble model consists of L bubbles connected 
in series, each bubble with w tubes. 

We choose the tube radii in the bubble model from the 
Hertz distribution 



/(r) = 2are~"'' 



(17) 



where a~^''^ is the characteristic pore radius. This form 
is often seen in experimental pore size measurements ||lq ] 
and has been used for modeling the pore size distribu- 
tion in filters p,pT|. For simplicity we assume identical 
tube lengths and measure downstream distance in units 
of the tube length, which is set equal to 1. We also 
assume that the flow rate in a tube of radius r is pro- 
portional to —r'^yP, where VP is the pressure gradient 
along the tube and fi depends on the nature of the flow, 
with /i = 4 corresponding to Poiseuille flow and /z = 2 to 
Euler fiow. Perfect mixing is assumed at each node. A 



A. Hydrodynamic dispersion in the bubble model 

In the large w limit, each bubble is nearly identical, and 
we can regard the particle motion as a directed random 
walk in which the average residence time r^. in bubble k 
{k = 0, 1, 2, . . .) is a random variable whose distribution 
R{t) is related to the flow induced entrance probability 
<I>(r) and the radius distribution /(r). This random walk 
description of the continuous particle motion introduces 
an additional stochasticity into the system. However, we 
will show below that this only modifies the hydrodynamic 
dispersion coefficient by an overall multiplicative factor. 

The master equations for pk{t), the probability that 
there is a particle in the fc"^ bubble at time t, are 

dpo _^ Po dpk Pk-i Pk ,^^s 



dt 



Po 

To' 



dt 



Tfc-l 



Pk_ 

Tk 



Here p is the initial particle number concentration, 
is the (constant) fiow rate, and the initial condition is 
Pfe(O) = for all k. Since the flow rate does not change 
in infiltration, constant pressure drop and constant fiow 
rate conditions are equivalent. 

The particle transport properties can be obtained in 
terms of the residence time distribution R{t), namely, the 
probability that a particle spends a time r in a bubble. 
This residence time distribution is related to microscopic 
distributions by 

R{r) = Hdr <i>{r) f{r) 6{t - -^), (20) 

Jo nr) 

where d{x) is the Dirac delta function. Since the fiow rate 
into a tube of radius r is 0$(r), the average flow velocity 
v{r) ~ 0$(r)/7rr^. Using this together with Eq. ( fil\ j for 
/(r), we obtain the flrst two moments of t 



(•00 

/ T'R{T')dT' = — 

Jo 'Pcf 

(r)2r(i + |)r(3-^) 



V 



(21) 
(22) 



where V = (wr^) is the average tube volume (recall that 
the tube length is fixed to be 1) and r(a;) is the gamma 
function. 

We solve Eq. ( [l9[) in the Appendix by the Laplace 
transform technique. From this solution, the average 
propagation velocity and the width of the front are 



Vo 



1 

Jt) 



2r(i + |)r(3-|) 



- 1 



(23) 



^{Dnt}^. (24) 



Thus we see that the dispersion coefficient is propor- 
tional to the average flow velocity (t)~^. When /i = 2 
(Euler flow) , the flow velocities in all the tubes are identi- 
cal and there should be no dispersion. However, Eq. ( |2^ ) 
gives a nonzero dispersion coefficient. As mentioned 
above, this arises from the stochasticity of the random 
walk picture for the particle motion. For the practically 
relevant case of /i ~ 4, the effect of this stochasticity is 
only to change the dispersion coefficient by a factor of 
order unity. 



B. Infiltration in the bubble model 

To describe infiltration in the bubble model, we need 
to specify the particle motion, the tube capacities, and 
particle trapping in a tube. For the particle motion we 
again assume that a particle chooses a tube according 
to flow induced probability and then advances with the 
average flow velocity v{r) of this tube. The capacity of 
a tube is proportional to its inner surface area, which 
is proportional to the tube radius, since all tubes have 
the same length. Last, the overall trapping probability 
of a tube is equal to the microscopic trapping probabil- 
ity 7 mult iplied by the fraction of open traps in a tube 
(Sec. |IIIB| ). 

To simulate this process efficiently we propagate the 
probability distribution function (PDF) of the suspended 
particles rather than simulating the motion of individual 
particles |19| , p0| . The PDF propagation therefore pro- 
vides the exact distribution of particle positions and tube 
capacities for a single realization of tube radii. Concep- 
tually, the PDF algorithm is equivalent to an exact inte- 
gration of the master equations. 

To implement the PDF propagation, we define 



that the average particle position is at the correct loca- 
tion along the tube, as illustrated in Fig. ^j. 

To construct the particle motion, let us temporar- 
ily disregard particle trapping. We set the time incre- 
ment to be St = 1/wmax, where Umax is the maximum 
fiow velocity among all tubes. In a time St^ a parti- 
cle at the entrance of the fastest tube should traverse 
the entire tube length which is equal to 1. We then let 
a particle in a slower tube, with velocity v^ < Umax, 



travel a distance 1 with probability u^ 



I 



main fixed with probability 1 



One can regard u^ 



as a normalized fiow velocity. By construction, such a 
particle travels the correct average distance in time 5t, 



i-< + o-(i 



<) 



Wfc/Wmax = Vl5t. 



Let us now recast this random walk into a probability 
propagation algorithm. Consider an element of the PDF 
which is at the junction before the /c*^ bubble. Before any 
particle motion occurs, we split this probability element 
among the downstream bonds in this bubble according 
to the flow induced probability at the tube entrance. We 
can view the probability element as advancing infinites- 
imally into each bond, as indicated on the left side of 
Fig. 0. Once this initial tube assignment is made, the 
probability element remains within its assigned tube un- 
til it reaches the next junction. 

Now consider the motion of a probability element 
which has just entered a particular bond. After a time 5t, 
a fraction ut of the PDF is advanced to the next bubble. 



while a fraction 1 
bond /3. 

t 



remains fixed at the entrance to 



t+dt 





„/3 



(2) (3) 



FIG. 5. Propagation of the PDF in tube j3 at bubble k. (1) 
Initial probability element p^{t). (2) Fraction remaining. (3) 
Fraction trapped. (4) Fraction advancing. This last element 
enters the next bubble and is then immediately split among 
the tubes according to the flow induced probabihties. 

Due to the filtration, a fraction of the flowing PDF 
becomes trapped in tube /3 in the fc"^ bubble at a rate 
which is proportional to the tube capacity c).{t). The 



pl{t): The probabilit,y that there is a particle at overafl trapping probabifity of this tube is therefore 

^fc = 7Cfe(0/Cfc(0)- After trapping has occurred, the 
tube capacity is decremented according to the following 
prescription. When one unit of PDF (equivalent to one 
particle) gets trapped, we dcflnc the bond capacity to 
be decreased by A. Therefore A is just the surface area 
covered by one particle. Correspondingly, c)^(0)/A equals 
the number of particles the tube can accommodate. 
Our algorithm for propagating an element of probabil- 



the entrance of tube /3 in bubble k (/3 = 1, 2, . . . , w, 
fc = 0,l,2,...). 

• c^(t): The capacity of tube (3 in bubble k. 

Since particles generally have different velocities, their 
positions could be anywhere within a tube. We simplify 
this by forcing particles to always be at the tube entrance 
by adjusting the time unit and the PDF propagation so 



ity at the entrance to bond [3 in the fc**^ bubble over a 
time 5t therefore consists of the following steps (Fig. g) : 

• Fraction of PDF remaining at the start: pf (1 — uf ). 

• Fraction trapped in tube (3: p^u^TJ^ . 

• Fraction advancing to the next junction: 

• Capacity change of the tube by trapping: 

The rate equations which account for these steps are: 
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(25) 
(26) 



The first term on the right hand side of Eq. (P5| ) is the 
fraction of probability that does not move, and the second 
term is contribution from elements of probability which 
has moved from the previous site. The flow induced prob- 
ability $f in Eq. (|5|) accounts for the fraction of PDF 
which enters into tube /3. 

To test this approach, we set 7 = (no trapping) in 
the above rate equations and simulate the PDF propa- 
gation. By this method, we find a traveling front whose 
basic properties coincide with the hydrodynamic disper- 
sion results given by Eqs. ( p3| ) and (|2J). 

It is also worth mentioning that our PDF algorithm 
can be generalized to allow for hopping a distance which 
is a fraction of the tube length. In this manner one can 
account for different longitudinal flow velocities at differ- 
ent radial positions within a tube (Taylor dispersion ||7|). 
In the limit of an infinitesimal hopping distance, continu- 
ous particle motion is reproduced by the PDF algorithm. 
Unfortunately, the gain in having a more accurate de- 
scription of the motion is offset by the complexity of the 
algorithm and large increase in the computation time. 



C. Asymptotic behavior 

To obtain the average invasion front profile over the 
tubes in each bubble, we define the bubble-average quan- 



We 



titles Pk{t) ^j^EppUt) and Qk{t) ^ ^Ep^it) 
first derive the invasion front velocity via the symmetry 
argument of Sec. |lIIB . A rigorous derivation of the front 
velocity from the master equations for Pk{t) and Qk(t) 
is given in |2l| ]. 

The carrier fluid is moving with velocity vq (Eq. (Eq)) 
and the input flux of invaders per tube is equal to VopV, 
since pV is equal to the number of invader particles per 
tube volume in the input fluid. On the other hand, in a 



reference frame moving with velocity vq, the input flux 
of defenders is equal to vqM where M = (r)/A is the 
average initial number of invaders a t ube ca n accommo- 
date. Following the argument in Sec. IIIB , we find the 
front velocities v and wo — w in the two reference frames 
are related by Mv — pV{vo — v), yielding 
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(27) 



Good filter performance means that the breakthrough 
time is long or, equivalently, that the propagation veloc- 
ity is slow. Eq. ( p7| ) implies that the propagation ve- 
locity can be made small by increasing the capacity of 
a pore, or by decreasing either the filter grain size or 
the input particle concentration. Notice that neither the 
reaction strength 7 nor the nature of the flow (through 
v{r) ~ r^~'^) affect this propagation velocity. 

We now study the asymptotic density profiles. Instead 
of working directly with the averaged quantities Pk{t) 
and Qk{t), we first focus on the behavior of a single tube 
of radius r, since tubes with the same radius in a bubble 
have identical time dependence. The asymptotic pro- 
files can be obtained after averaging over the distribu- 
tion of tube radii. Therefore, we label tubes according 
to their radii instead of the index (3. We denote Pk{t; r) 
and Ck{t; r) as the PDF and capacity of a tube of radius 
r in the fc*^ bubble. 

Let us first focus on the PDF profile in the invaded 
region. Here, traps are mostly unoccupied, so that the 
tube capacity Ck{t; r) is approximately equal to its initial 
value Cfc(0;r), which is proportional to the tube radius 
and we set it equal to r. The arbitrariness in the unit 
of capacity can be controlled by the magnitude of the 
parameter A. Then Eq. (|25|) becomes. 



Pk{t + St; r) ~ pk(t] r) (1 - u{r)) 

+ (l-7)$(r) fdr'fir')pk-i{t;r')uir'), (28) 



where the integration over r' replaces the summation over 
the tube index. The flow induced probability <I>(r) and 
normalized velocity u{r) are independent of the down- 
stream position k because these only depend on the ra- 
dius of a tube. 

The equation for P{x, t) can be obtained by multiply- 
ing Eq. ( pq ) by /(r) and integrating over r. As in Sec. Ill, 
we take 7^1 and consider the continuum limit. If we 
redefine the length of the bubble from 1 to Sx, St becomes 
Sx/vjnax- Integrating Eq. ( p8|) and expanding in St and 
Sx yields 



StdtP{x, t) = - I drf{r)u{r)[Sxdxp{x, t; r) + jp{x, t; r)] 

(29) 

Here, we use Jdrf{r)^{r) — 1. Dividing Eq. ( p9| ) by St 
changes u{r) back to v{r). After redefining 7 — > 7/fe as 
before, we obtain 



dtP^- drf{r)v{r) [OM^, *; + IPi^, t; r)] . (30) 



In the steady state co-moving frame, Eq. (pO^ becomes 



•jdeP 



drf{r)v{r) [d^p{^; r) + 7p(^; r)] 



(31) 



Since only a small number of particles have entered 
the invaded region, the density of moving particles is 
approximately proportional to r**, and we introduce the 
ansatz p(^;r) = r'^g{£^) to factorize the PDF. In order to 
calculate the dominant contribution from the integral in 
Eq. (31^), we substitute v{r) =110-1- Sv{r), where Sv{r) is 
the deviation from the average carrier fluid velocity vq. 
Since Sv{r) has zero mean, the dominant contribution to 
the integral over v{r) in Eq. (31) comes from the con- 
stant part vo- Using these approximations in Eq. (31), 
and using P{^) = (?'^).9(0; "^^ fi^^d 



vg' ~vo(g' + jg). 



(32) 



Since v < vo, we find g(^) ~ exp [—vqj^/{vo — v)]. Hence 
the profile of free particles in the invaded region P{£,) de- 
cays exponentially in ^, with a characteristic decay length 
which has the same I/7 dependence as in the ID model. 
Let us now turn to the analysis of the capacity pro- 
file in the tail region. In terms of p(a;,t;r) and c{x,t;r), 
Eq. ( p6| ) becomes 

dtc{x,t;r) = ^Ajpix,t;r)v{r)^^^^^. (33) 
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where the last relation serves to define the profile expo- 
nent X. This is one of our primary results. Correspond- 
ingly, the PDF in the tail region will approach its initial 
value with the same power law. 

The existence of the power law tail in the capacity pro- 
file stems from the fact that the flow rate is not affected 
by trapping. Thus when large pores are "used up" , the 
fluid still predominantly flows through these pores, lead- 
ing to a substantial unused capacity in the smaller tubes. 
It is these unused smaller tubes which contribute sub- 
stantially to the capacity profile in the tail region. This 
mechanism is quite general and only depends weakly on 
the form of the radius distribution. For example, for a 
uniform distribution in the range r — (0, 1), we obtain 
A = 2/(/i— 1). However, if there is a finite lower cutoff in 
the radius distribution, the PDF will have an asymptotic 
exponential tail. 

It is interesting to note that the density profile has dif- 
ferent dependence on 7 in the invaded and tail regions. 
From Eq. (^) , the density profile contains an overall fac- 
tor 7^^. Thus 7 typically does not appear as an overall 
scale factor of the entire profile, as in the invaded region. 
However, for the practically relevant case of fi — 4, the 
exponent in Eq. ( |37| ) is equal to 1, and 7 becomes the 
overall scale factor of the profile. 



Since there is negligible trapping in the tail region, the 
particle motion follows that of the carrier fluid. Thus 



p{^;r) 



Trr^p, 



where 



is the tube volume, and the 



flow velocity is v{r) = (/)4>(r)/7rr^ = (f)r^~'^ / -k {r^) . Sub- 
stituting these in Eq. ( |33| ) and transforming into the co- 
moving frame gives 



d^c[S,]r) ~ sr'' ^c[S,]r), 



(34) 



where s = A7p(/)/i;(r^) denotes the strength of the parti- 
cle trapping reaction. We now integrate Eq. (pj) from 
—^1 <C to ^2 ~ and use the boundary condition 
c{^2, r) ~ r to obtain 
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where we drop the subscript of ^1 . 

Finally, the average bond capacity as a function of po- 
sition with respect to the front, Q{C), is 
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(36) 



For large |^|, the integral is dominated by the smallest 
tubes and the initial distribution of tube radii is irrel- 
evant in the tail region. Hence the factor ar^ in the 
exponential can be ignored. Performing the resulting in- 
tegration gives 



D. Numerical results 



In our numerical simulations, we set the input particle 
flux per tube pcj) equal to 1, which means that w units 
of PDF are injected into the system at every time step. 
This can be achieved by choosing p ~ I/tt and cj) — n, 
which also makes vq = a (Eqs. ( pi|) and (p3|)). 
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FIG. 6. Normalized profiles for a single realization. Pa- 
rameters used are p. — 4, a = 1, 'y = 0.1, and A — 0.4. Gray 
lines are raw data, black lines are smoothed. 



We applied the PDF propagation of Eqs. ( |25| ) and (|2£ 
to a system of size w x L = 200 x 1024. Due to the exact 
nature of the PDF algorithm, a single realization provides 
good quality data for w = 200 tubes. A system length 
of L = 1024 is sufRciently long to give the continuum 
fmictional form of the profiles. All the data shown be- 
low are results of single realization of tubes including the 
network simulation in the next section. The simulation 
is stopped before the front exits the system. 

Density Profiles. Fig. ^ shows typical particle and tube 
capacity profiles. There are strong bubble-to-bubble fluc- 
tuations and some type of smoothing procedure is nec- 
essary. We use the Savitzky-Golay smoothing technique, 
which approximates successive windows of data points to 
a 4'^ order polynomial (solid lines in the figure) ||2^ . This 
technique is superior to local averaging because Savitzky- 
Golay smoothing can faithfully follow rapid changes in 
the profile, as can be seen in Fig. ||. This smoothing is 
also useful in estimating the exponents. Since logarithm 
of the profile in the tail region amplifies the fluctuation in 
nonlinear way, slopes of the raw data in Fig. are larger 
than those from the smoothed data, and differs from the 
predicted value of the profile exponent A. 
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FIG. 7. Front position vs. time for different microscopic 
parameters for (a) the bubble model and (b) the square lat- 
tice. Filled circles are for the same parameters as those in 
Fig. y. Other data modify these parameters as indicated in 
the legend. The straight lines are linear fits to the data with 
slopes as (a) 0.33, 0.32, 0.31, 0.19, 0.12, and (b) 0.32, 0.32, 
0.31, 0.19, 0.12, from top to bottom. The corresponding ve- 
locities predicted by Eq. (||) are 0.311, 0.311, 0.311, 0.184, 
and 0.119. 

Front Velocity. Fig. R(a) shows the front position, de- 
fined as the point where Pk{t) is half of its saturation 
value, versus time. Notice that a constant front propaga- 
tion velocity sets in almost immediately. With p = 1/tt 
and 4> — 1^, Eq. (p7|) gives 



1 



vq 



(38) 



No- 



+ V7ra/2A 

The slopes in Fig. ^(a) agree well with Eq. ( ^ 
tice that the propagation velocity does not depend on 
the reaction strength 7 nor the exponent /i in the radius 
dependence of the velocity. 



Tail Profile. Fig. ||(a) shows the tube capacity profile 
Q{i) in the tail region as a function of the distance |^| 
(^ < 0) from the front on a double logarithmic scale. The 
plot becomes straight for large |^| and the slope in this 
region corresponds to the exponent A = 3/(/i — 1) pre- 
dicted by Eq. (|37|). For the uniform distribution on (0, 1), 
we predicted the profile exponent to be A = 2/(/i — 1). 
For /i = 4, the exponent value of 2/3 agrees well with our 
simulations (Fig. 0). However, for a radius distribution 
with a lower size cutoff, we expect an exponential density 
profile (inset to Fig. pi). 
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FIG. 8. Tube capacity profile in the tail region on a double 
logarithmic scale for /i = 4, 3, and 2 (top to bottom). Other 
parameters are the same as in Fig. |6[ (a) Bubble model: the 
slopes of the data in straight region are 0.97 (1), 1.49 (1.5), 
and 3.17 (3). The numbers in parenthesis are the prediction 
3/(/i — 1). (b) Square lattice: Solid lines are smoothed data. 
The thick straight lines are linear fits, with slopes 0.95, 1.07, 
and 1.28. 




FIG. 9. Tube capacity profile on a double logarithmic scale 
for a uniform radius distribution on (a, 6). Parameters other 
than the radius distribution are the same as in Fig. 0. The 
slope of the straight region for (a, b) — (0, 1) is 0.63 (predicted 
value 2/3). Inset: log-normal plot when (a, b) = (0.3, 1.3). 



As we also discussed in Sec. [V C, the amplitude of the 



density profile in the tail region typically has a power- 



law dependence on 7. For a Hertz distribution of parti- 
cle radii, this amplitude should be proportional to 7^ 
according to Eq. @. Thus Fig. |l^(a) shows QiO^^^l^l 
versus l/\£,\ for /i = 4 and /i = 3. Values of the abscissa 
should be proportional to I/7, which is indeed the case. 
Invader Profile. Fig. |ll|(a) is a log-normal plot of the 
invader profile P{^) (raw data) versus ^. The slopes of 
the two data sets are in excellent agreement with the pre- 
dicted value wo7/(wo — v) from Eq. (p^). In the invaded 
region, since there are few particles present, the corre- 
sponding PDF monotonically decreases and its fluctua- 
tion is significantly smaller compared to the tail region. 
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FIG. 10. Dependence of the capacity profile in the tail re- 
gion on 7. Only smoothed data are shown. Solid curves: 
fj, = 4, dashed curves: ^ = 3. Thick curves: 7 — 0.05, thin 
curves: 7 = 0.1. Other parameters are the same as in Fig. pi 
(a) Bubble model: The value A = 3/(/x — 1) is used, (b) Square 
lattice: From Fig. ||(b), A = 0.95 for /i = 4 and 1.07 for ^ = 3 
are used. Notice that the values of abscissa are proportional 
to 1/7 from Eq. (pTI). Arrows are guides to the eye. 
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FIG. 11. Log-normal plot of the invaded region. The same 
data sets as in Fig. flG are used. The slopes of the straight 
regions are, from left to right, (a) bubble model: 0.152, 
0.150, 0.0744, 0.0748, and (b) square lattice: 0.152, 0.155, 
0.0756, 0.0773. The corresponding values of the slope from 
uo7/(uo - v) are 0.145 and 0.0726. 



V. SQUARE LATTICE NETWORK 

We now consider infiltration on a square lattice net- 
work of tubes. Here, local mixing at tube junctions oc- 



curs as opposed to the mean-field-like mixing in the bub- 
ble model. Nevertheless, many of our predictions from 
the bubble model continue to be valid for the lattice net- 
work. For example, we expect that the propagation ve- 
locity given by Eq. ( p7| ) will continue to hold in the lat- 
tice network because it is determined only by the bound- 
ary conditions of the particle flux and initial filter capac- 
ity. We also find numerically that the network model has 
both the same exponential invader profile and the power 
law capacity profile as in the bubble model, although the 
values of the amplitudes and decay exponents are differ- 
ent. Overall, it appears that the bubble model provides 
an excellent account of the numerical results from the 
lattice network. 



A. The model 

We study a square lattice of size w x L — 200 x 1024 
which is tilted at 45°. A periodic boundary condition is 
imposed in the transverse direction. The tube radii are 
drawn from the Hertz distribution of Eq. (J17|). Notice 
that the bubble model would arise from this system by 
merging together all sites at the same longitudinal posi- 
tion. The overall fiow rate cj) is set to tt and the particle 
density to p = I/tt, just as in the bubble model. 

For a given set of tube radii, the flow field is calculated 
by using the conjugate gradient method ||2^ to solve the 
set of linear algebraic equations for fluid conservation at 
each node. The tolerance of the computation is set so 
that the measured average PDF and the tube capacity 
in the fc"^ layer are accurate to within 0.01%. After the 
flow fleld is solved, we use the same PDF algorithm as in 
the bubble model to track the motion of the suspended 
particles. The new features due to the lattice nature of 
the network are that tubes are only locally connected and 
that the local flow direction is not always downstream. 



B. Numerical results 

To facilitate comparisons, all our numerical results for 
the bubble model and the square lattice are presented 
side-by-side. Fig. u\ shows the front position versus time 
for both the bubble model and the square lattice. The 
square lattice results are in excellent agreement with the 
bubble model prediction for the front velocity, Eq. Pq). 
Similarly, the tube capacity profile exhibits a power law 
tail (Fig. |g(b)). However, the dependence of the decay 
exponent A on ^ is much weaker compared to the bubble 
model. As /z decreases from 4 to 2, A increases slowly 
from 0.95 to 1.28. 

To isolate the dependence of the capacity profile on the 
trapping probabihty 7, Fig. |l^(b) shows Q{S,)^^^\^\ ver- 
sus 1/ICI in the tail region. Here, the values of A used are 
obtained from Fig. @(b). Unlike the bubble model, the 
overall amplitude has a relatively stronger dependence 
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on fx. However, the dependence on I/7 still holds even 
in the network case. It seems that in the square lattice 
network, /j, has more effect on the amplitude of the tail 
than on the decay exponent. 

Lastly, Fig. l^(b) shows the density profile of invaders 
in the invaded region. The slopes of these particles are 
almost identical to those of the bubble model. Currently 
we do not have clear explanation of this fact. It seems 
that the network geometry does not affect the profile. In 
fact, the characteristic decay length, (vq — u)/wo7, of the 



profiles in Fig. 1 1 is the same as that of the correspond- 



ing ID model in Sec. [II B where pV invaders are injected 
with velocity vq into the chain of defenders of capacity 
M. 



VI. SUMMARY AND DISCUSSION 
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FIG. 12. Qualitative dependence of the output concentra- 
tion curve (a) on the propagation velocity, and (b) on the 
reaction strength. 



In this paper, we studied infiltration, in which sus- 
pended particles are removed from a carrier fluid as the 
suspension passes through a porous medium. The trap- 
ping mechanism has a built-in saturation so that once all 
available trapping sites are used up, subsequent particles 
can pass through the medium freely. The particles are 
assumed to be sufficiently small that their trapping does 
not change the flow rate. The basic dynamical proper- 
ties of this infiltration process are the density profile of 
the invader particles and the capacity profile of the re- 
maining active pores. When the invader profile reaches 
the end of the system, the output concentration of parti- 
cles quickly increases to a saturation level and the filter 
should be discarded. Thus the features of this profile are 
important to understand the operating characteristics of 
infiltration. 

We have developed a series of discrete network mod- 
els to describe the basic characteristics of infiltration, 
starting with a one-dimensional model and building up 
to the bubble model, which is a series array of parallel, 
multiple-capacity tubes. The advantage of these quasi- 
onc-dimcnsional models is that they remain relatively 
simple, even after incorporating local spatial heterogene- 
ity. The bubble model, in particular, appears to capture 
many of the quantitative features that we observed in nu- 
merical simulations of infiltration on a square lattice tube 
network. Our modeling is also fiexible, so that variations 
can be easily implemented for case-specific situations. 

Our main qualitative result is that basic dynamical 
features of the system, including the value of the front 
propagation velocity, the exponential profile of flowing 
particles in the invaded region, and the power law ca- 
pacity profile of pores in the tail region, are relatively 
insensitive to microscopic details of the model. We have 
also identified the basic parameters which do affect quan- 
titative features of the profiles. It is useful to summarize 
these results and to compare with experimental data, as 
well as with predictions from previous studies based on 
the reaction-diffusion equation approach |0,@,E3l. 



Invader Concentration at the Output. In typical exper- 
iments, the invader concentration at the output is mea- 
sured as a function of time. A slower propagation velocity 
shifts this output concentration curve to a later time, as 
indicated in Fig. [L2|(a). Typically, the time unit is nor- 
malized by the time for passive particles to pass through 
the system [||Q] . Hence, the amount of the time shift is 
determined by the ratio v/vq rather than by v itself. 

A nice set of infiltration experiments, as well as an 
accompanying numerical study of the reaction-diffusion 
equation were performed in [pi. In these experiments, 
contaminant solutions with different values of the invader 
particle diameter d, but with fixed mass concentration 
were used. This makes the corresponding number density 
p of invaders in each solution proportional to d~^. Also, 
since the cross-sectional area of each particle is propor- 
tional to d^ , the average initial tube capacity M varies as 
d"2_ FromEq. (|2^), we then have w/wo -- (H-d)"^ Thus 
the output concentration curve shifts to a later time for 
a solution with a larger value of d, which is consistent 
with the experiment and numerical predictions in ||6|]. 

In another set of experiments, different electrolyte con- 
centrations of the carrier fluid were used. This mainly 
affects the trapping rate 7. For a larger trapping rate, 
the width of the output concentration curve, namely, the 
time range over which the output concentration changes 
from zero to its saturation value becomes narrower. How- 
ever, there is no shift in the breakthrough time because 
the propagation velocity is independent of 7. These two 
features are illustrated in Fig. |2|(b). This behavior again 
qualitatively agrees with the experiment in [pj . 

It would also be interesting to study the effect of dif- 
ferent fllter grain sizes. In our model this would be ac- 
complished by changing the characteristic parameter a 
of the pore size distribution. In turn, this affects the in- 
vader propagation velocity, and the output concentration 
curve will shift in time accordingly. However, since the 
microscopic parameters we use in our modeling may be 
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coupled with each other in experimental situations, dif- 
ferent sizes of filter grains or invader particles may also 
affect, e. g., the reaction strength 7. We can incorporate 
such a coupling effect by extending our model to deal 
with these effects explicitly. For example, we can adapt 
the microscopic models of particle trapping on a single 
sphere or on a plane [§J8), to the tube geometry. From 
such an approach, we can express the reaction strength 
7 as a function of the invader or defender diameter. 

Tail of the Output Concentration Curve. A slowly de- 
caying tail in the deviation of the output concentration 
from its asymptotic value is generally observed in exper- 
iments |0,p|,p3|. This observation is in contrast to the 
empirical approaches, such as that given in M, which 
gives an exponential profile for the whole time range. 
Their prediction agrees with experimental observations 
at early times but then deviates at later times, implying 
that the output profile at later times is not exponential. 

A closely related approach, based on the study of a 
reaction-diffusion equation, is presented in Q , along with 
experiments which measure the output concentration. 
Unlike Q, here the adsorption rate depends on the lo- 
cal concentration of contaminants and thus is spatially 
inhomogeneous. A crossover from a rapid increase to a 
slowly decaying tail of the output concentration was nu- 
merically predicted. However, the functional forms of 
these two regimes - in particular, whether they are ex- 
ponential or power law in time - were not investigated 
quantitatively. However, the data presented in this work 
seem consistent with a slower than exponential decay of 
the density profile. 

In [p3| , an exponential output concentration profile, 
c{x) ~ exp(— Ax), is assumed from the outset, where 
X is the downstream distance and A is the experimen- 
tally measured filter coefficient. The corresponding ex- 
perimental data show that A is constant at early times, 
and then sharply decreases at later times. Thus the ini- 
tial stage of the experiment is consistent with an expo- 
nential profile, but later the profile decays more slowly. 
In |2^ , this is attributed to a "blocking effect" in which 
previously-deposited particles can block the further de- 
position of particles onto nearby available trapping sites. 

Probability of Encountering an Open Trap. As a last 
remark, let us examine the assumption that the prob- 
ability of encountering an open trap is proportional to 
the fraction of open traps in a pore (Sec. [II B| ). Sup- 
pose instead that one takes into account the volumetric 
effect that particles far from the surface of the pore do 
not have chance to encounter a trap. Then the fraction 
of particles in contact with the inner surface of a tube is 
proportional to 1/r, namely, the ratio between the sur- 
face area and the volume of the tube. This would lead 
to the interaction terms involving T^ in Eqs. (ESh and 
( Pq ) being multiplied by another factor of 1/r. However, 
this modification does not affect the propagation veloc- 
ity of the front, nor the power law feature of the tail. 
Only the decay exponent changes through the steps of 
Eqs. (p3)-(pfl) with an additional factor of 1/r. 



Our results can also provide practical guidelines for 
improving the design of a filter in two aspects, namely, 
the breakthrough time and the amount of filter material 
used before the breakthrough. A longer breakthrough 
time can be achieved by having a smaller filter grain size, 
a lower input concentration, or a larger pore capacity. 
While these trends may seem intuitively clear, we can 
quantitatively estimate the increase in the breakthrough 
time thr oug h the expression for the propagation veloc- 
ity, Eq. (p7|). When breakthrough occurs, the amount of 
unused filter material is determined by the shape of the 
tail in the density profile of the defenders. According to 
Eq. (^), the amplitude of this tail is proportional to 7"'*'. 
From this, we can quantitatively estimate the amount of 
filter material left unused at the breakthrough time as a 
function of the reaction strength. 
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{p^) = l^{l-B)B\ B^ 



1 



(A5) 



where (•) denotes an average over R{t). Note that B 
does not depend on m because all the r^'s are indepen- 
dent and identically distributed in the large w limit. For 
the long-time limit, we need the small-s behavior of B. 
Accordingly, we expand B = 1 — s{t) + s^ {t^) ... in terms 
of the moments (r™). 

The profile of the carrier fluid (dashed line in Fig. g) 
is monotonically decreasing near the front, thus can be 
characterized by its first derivative, which gives a bell 
shaped distribution centered at the front. We divide the 
derivative by the total sum of derivatives X]m((P™) ~ 
(Pm-i)) = ~(po) to obtain the normalized probability 
distribution of the front 



APPENDIX A: SOLUTION OF THE 
HYDRODYNAMIC DISPERSION EQUATION 

We solve the master equations, Eq. (|l9), by 
the Laplace transform method. Define Pk{s) = 
/p°° dtpk{t)e~'^*, and take Laplace transform of Eq. (|l9) 
to find 



spo 



Po 



S To 

Rearranging yields 
Po, 



spk 



Pk-i _ Pfc 

Tfe-l Tk 



To 



-(1 + ros) 



Pfe.-, ^ . Pk-i 

(1 +TkS) = 

Tk Tk-1 



(Al) 



(A2) 



V{k) 



1 



{Po) 



[{Pk- 



{Pk 



(A6) 



The average position of the front is, using Eq. (|A5| ), 
and the steady state solution of Eq. (p^), (po) = p(t){T), 

l-B 



^^HV^EMi?'-^-^'] 



(A7) 



We multiply the above equations for indices 0, 1, 2, • • • k 
together and rearrange terms yet again to obtain 



(1 +TmS) = . 

m—0 



(A3) 



Before taking averages, we use Tk = ^ [1 — 1/(1 + r^s)] (1- 
Tks) on the left hand side of Eq. (A3) to get 



where the over-bar means averaging over 'P{k), and C ^ 
is the inverse Laplace transform. We use the identity 
£-i(s-/5) = tl^-'^/T{l3) for the last step. We assume suf- 
ficiently long chain of bubbles in the summations above 
to prevent finite length effect. From above, we find the 
propagation velocity as w = (r)^^ (Eq. (p3|)). 

Similarly, (P) = (po)"^ Y.k ^^[(Pfe-i) - {Pk)], and the 
width of the front is 



PkJ{{l 

m=0 



TmS) = ^r 1 - 



1 



l+TfeS 



(l + Tfes). (A4) 



Averaging over the residence time distribution R{t) de- 
fined by Eq. (EQ), we obtain 



2, 1 



[{k^)^{k)y^i — 



t 



2r(i + ^)r(3-^) 



(A8) 



which gives Eq. (P4| ) 
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